#############################################################################################################################################################################################################
### Setup

# Set working directory
setwd('~/Dropbox/Projects/Minority Party/Paper 2/APSR Replication Files/')

# Source setup file
source('setup.R')


### Load bills data
bills <- fread(file='bills97-109_partyPriorities_issuePreferences.csv', data.table=FALSE, stringsAsFactors = FALSE) %>%
  filter(congress %in% 99:109)


### Add in SOTU data
# Load, filter to years matching the congresses studied (1985-2006---though the 109th ends in 2007 it was before the 2007 sotu)
sotu <- fread(file='SOTU-topicmentions-1946-2019.csv', data.table=FALSE, stringsAsFactors=FALSE) %>%
  filter(year %in% 1985:2006) %>%
  dplyr::rename(minor=subtopic) 

# Match years to congresses
sotu$congress <- NA
congs <- sort(unique(bills$congress))
years <- sort(unique(sotu$year))
for(ii in 1:d(years)){
 if(years[ii] %% 2 == 1){ # If it's an odd year
  sotu$congress[sotu$year==years[ii]] <- congs[(ii+1)/2]
 } else {
  sotu$congress[sotu$year==years[ii]] <- congs[(ii/2)]
 }
}

# Drop years
sotu %<>% dplyr::select(-year)

# Join sotu with bills
bills <- left_join(bills, sotu, by=c('congress', 'minor'))

# Add in 0s for sotu variables (issues that weren't mentioned)
bills$any_statement[is.na(bills$any_statement) & !is.na(bills$minor)] <- 0
bills$num_statements[is.na(bills$num_statements) & !is.na(bills$minor)] <- 0
bills$perc_statements[is.na(bills$perc_statements) & !is.na(bills$minor)] <- 0
bills$perc_topic_statements[is.na(bills$perc_topic_statements) & !is.na(bills$minor)] <- 0

# Cut down to distinct rows (if any duplicates)
bills %<>% distinct()

# Indicator for divided government
bills %<>% mutate(divided=ifelse(!congress %in% c(103, 107, 108, 109), 1, 0))

# Remake majPriorityAll and minPriorityAll into factors (makes plotting better)
bills$majPriorityAll <- factor(bills$majPriorityAll, labels=c('No', 'Yes')) %>% relevel(., ref='No')
bills$minPriorityAll <- factor(bills$minPriorityAll, labels=c('No', 'Yes')) %>% relevel(., ref='No')
bills$any_statement <- factor(bills$any_statement, labels=c('No', 'Yes')) %>% relevel(., ref='No')

# Filter to substantive bills in the 99th-109th congresses
modDat <- filter(bills, substantive==1)

# Define plot colors
cols3 <- brewer.pal(6, 'Greys') %>% .[c(4:6)]



#########################################################################################################################################################################
### Regression models

# Full model (opportunity and cohesion), no interaction
m.noint <- glm(pv_all ~ majSD + minSD + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll +
                sizeMajH + numSenate + unityPty + presControl, 
              family=binomial, data=modDat)

# Full model (add interaction)
m.twoint <- glm(pv_all ~ majSD * minSD + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll + 
                sizeMajH + numSenate + unityPty + presControl, 
              family=binomial, data=modDat)

# Model including SOTU mentions--adds motivation analysis 
# (only in cases of divided government between House and White House, which cuts out 103rd, 107th-109th congresses)
m.sotu <- glm(pv_all ~ majSD * minSD * any_statement + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll +
                sizeMajH + numSenate + unityPty, 
              family=binomial, 
              data=filter(modDat, divided==1))

# Models computing whether a bill became law (among all bills and those that passed the House)
covs <- ' + majority + leader + les + mref + absPVI +
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll +
                sizeMajH + numSenate + unityPty + presControl'
m.law <- glm(paste0('plaw ~ majSD * minSD * any_statement', covs),
             family=binomial, data=filter(modDat, divided==1))
m.law.pass <- glm(paste0('plaw ~ majSD * minSD * any_statement', covs),
                  family=binomial, data=filter(modDat, divided==1 & passh==1))

# Compute R^2 for those going in the text table
r2s <- sapply(list(m.noint, m.twoint, m.sotu, m.law, m.law.pass), 
              function(x){
                round(NagelkerkeR2(x)$R2, 3)
              })



################################################################################################################################################################################################################
### Table 1: Regression models

# Save model output (slight differences in formatting between this and the paper due to availability in stargazer options)
stargazer(m.noint, m.twoint, m.sotu, m.law, m.law.pass, type='latex', out='table1.tex',
          column.sep.width='0pt', no.space=TRUE, font.size='tiny', dep.var.labels.include = FALSE, dep.var.caption='', label='regTab',
          title='Logistic regression models predicting whether a bill receives a final passage House vote (columns 1-3) or whether a bill becomes law (columns 4 and 5).', model.numbers=FALSE,
          column.labels = c('Agenda (1)', 'Agenda (2)', 'Agenda (3)', 'Lawmaking (4)', 'Lawmaking (5)'),
          covariate.labels = c( # Main IV labels
                               'Majority Spread', 'Minority Spread', 'Minority Priority',
                               # Sponsor covariate labels
                               'Majority', 'Leadership', 'Legislative Effectiveness', 'Member Referral Committee', 'Electoral Safety', 
                               # Bill covariate labels
                               'Closer to Minority', 'Election Year', 'Related in Senate', 'Party Difference', 'Majority Priority', 
                               # Political environment covariate labels
                               'House Majority Size', 'Seats in Senate', 'Party Unity', 'White House Control', 
                               # Interaction labels
                               'Majority Spread * Minority Spread', 'Majority Spread * Minority Priority', 'Minority Spread * Minority Priority', 
                               'Majority Spread * Minority Spread * Minority Priority',
                               # Constant
                               'Constant'
                               ),
          keep.stat='n', add.lines=list(c('Nagelkerke Pseudo R2', r2s))
          )



################################################################################################################################################################################################################
### Figure 3: 2-way interaction (majority spread and minority spread)

# Get data
pDat <- predDat(model=m.twoint, vars=c('minSD [all]', 'majSD'), data=modDat, std=TRUE)

# Plot
pdf("figure3.pdf", width=5, height=5)

breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.08), xaxt='n', las=2, main='')
for(ii in 1:nrow(pDat)){
  colIndex <- which(c('Low', 'Median', 'High')==pDat$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pDat$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pDat$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pDat$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

# Turn off device
dev.off()



################################################################################################################################################################################################################
### Figure 4: 3-way interaction (majority spread, minority spread, minority priority)

# Get data
p3 <- predDat(model=m.sotu, vars=c('any_statement', 'minSD', 'majSD'))

# Set up plot parameters
pdf("figure4.pdf", width=10, height=5)
par(mfrow=c(1,2))

# Plot
p3No <- filter(p3, x==1)
breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.08), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(p3No)){
  colIndex <- which(c('Low', 'Median', 'High')==p3No$majSD[ii])
  if(ii %in% 1:3){
    points(ii, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

p3Yes <- filter(p3, x==2)
breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.08), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(p3Yes)){
  colIndex <- which(c('Low', 'Median', 'High')==p3Yes$majSD[ii])
  if(ii %in% 1:3){
    points(ii, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topleft', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

# Reset, turn off device
par(mfrow=c(1,1))
dev.off()



################################################################################################################################################################################################################
### Figure 5: Probability of becoming law, 3-way interaction 

# Get data
pLaw <- predDat(m.law, vars=c('any_statement', 'minSD', 'majSD'))

# Set up plot parameters
pdf("figure5.pdf", width=10, height=5)
par(mfrow=c(1,2))

# Plot
pLawNo <- filter(pLaw, x==1)
breaks <- 3
plot(0:(nrow(pLawNo)+2*breaks), 0:(nrow(pLawNo)+(2*breaks)), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0,0.04), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(pLawNo)){
  colIndex <- which(c('Low', 'Median', 'High')==pLawNo$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pLawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pLawNo$lo[ii], pLawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pLawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pLawNo$lo[ii], pLawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pLawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pLawNo$lo[ii], pLawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

pLawYes <- filter(pLaw, x==2)
breaks <- 3
plot(0:(nrow(pLawYes)+2*breaks), 0:(nrow(pLawYes)+(2*breaks)), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0,0.04), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(pLawYes)){
  colIndex <- which(c('Low', 'Median', 'High')==pLawYes$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pLawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pLawYes$lo[ii], pLawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pLawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pLawYes$lo[ii], pLawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pLawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pLawYes$lo[ii], pLawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

# Reset, turn off device
par(mfrow=c(1,1))
dev.off()



################################################################################################################################################################################################################
### Figure 6: Probability of becoming law (among bills passing the House), 3-way interaction

# Get data
pLawPass <- predDat(m.law.pass, vars=c('any_statement', 'minSD', 'majSD'))

# Set up plot parameters
pdf("figure6.pdf", width=10, height=5)
par(mfrow=c(1,2))

# Plot
pLawPassNo <- filter(pLawPass, x==1)
breaks <- 3
plot(0:(nrow(pLawPassNo)+2*breaks), 0:(nrow(pLawPassNo)+(2*breaks)), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0.2,0.6), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(pLawPassNo)){
  colIndex <- which(c('Low', 'Median', 'High')==pLawPassNo$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pLawPassNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pLawPassNo$lo[ii], pLawPassNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pLawPassNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pLawPassNo$lo[ii], pLawPassNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pLawPassNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pLawPassNo$lo[ii], pLawPassNo$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

pLawPassYes <- filter(pLawPass, x==2)
breaks <- 3
plot(0:(nrow(pLawPassYes)+2*breaks), 0:(nrow(pLawPassYes)+(2*breaks)), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0.2,0.6), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(pLawPassYes)){
  colIndex <- which(c('Low', 'Median', 'High')==pLawPassYes$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pLawPassYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pLawPassYes$lo[ii], pLawPassYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pLawPassYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pLawPassYes$lo[ii], pLawPassYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pLawPassYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pLawPassYes$lo[ii], pLawPassYes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='bottomright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

# Reset, turn off device
par(mfrow=c(1,1))
dev.off()



######################################################################################################################################################################
### Use passage as the dv instead of reaching a vote

### Set up data

# Filter to substantive bills in the 99th-109th congresses
mdat <- filter(bills, (congress %in% 99:109) & substantive==1)

# Vectors of IVs and covariates
covs <- c('minPct', 'closerMinCS', 'seniority', 'leader', 'les', 'mref', 'elec.year', 'relSen', 'sizeMajH', 
          'numSenate', 'unityPty', 'senControl', 'presControl', 'absPVI') 
ivPref <- c('majSD', 'minSD', 'medDiff') 
ivPri <- c('majPriorityAll', 'minPriorityAll')


### Run final models

# Model formulas
modForm <- formula(paste0('passh ~ ', paste(ivPref, collapse='*'), ' + ', paste(ivPri, collapse='*'), ' + majority + ', paste(covs, collapse='+')) )
modFormSub <- formula(paste0('passh ~ ', paste(ivPref, collapse='*'), ' + ', paste(ivPri, collapse='*'), ' + ', paste(covs, collapse='+')))
modFormNoPriInt <- formula(paste0('passh ~ ', paste(ivPref, collapse='*'), ' + ', paste(ivPri, collapse='+'), ' + ', paste(covs, collapse='+')))
modFormInt <- formula(paste0('passh ~ ', paste(ivPri, collapse='*'), ' * ', paste(ivPri, collapse='*'), ' + majority + ', paste(covs, collapse='+')))

# Run models
m.twoint <- glm(modForm, family=binomial, data=mdat)
m.noPriInt <- glm(modFormNoPriInt, family=binomial, data=mdat)
m.min <- glm(modFormSub, family=binomial, data=filter(mdat, majority==0))
m.maj <- glm(modFormSub, family=binomial, data=filter(mdat, majority==1))
m.int <- glm(modFormInt, family=binomial, data=mdat)

# Summaries/fit
summary(m.twoint)
summary(m.min)
summary(m.maj)
summary(m.int)
NagelkerkeR2(m.twoint)
NagelkerkeR2(m.min)
NagelkerkeR2(m.maj)
NagelkerkeR2(m.int)

# Compute R^2 for those going in the text table
r2.full <- NagelkerkeR2(m.twoint)$R2 %>% round(., 3)
r2.min <- NagelkerkeR2(m.min)$R2 %>% round(., 3)
r2.maj <- NagelkerkeR2(m.maj)$R2 %>% round(., 3)

# Save model output
stargazer(m.twoint, m.min, m.maj, type='latex', 
          out='~/Dropbox/Projects/Minority Party/Analysis/paper2_regTable_appendix.tex')


### Predicted Probabilities

# Write a function
predDat <- function(model, vars, sd=FALSE){
  # Get the predicted data from plot_model
  p <- plot_model(model, type='pred', terms=vars)
  
  # Pull out the predicted data
  dat <- as.data.frame(p$data)
  
  # Pull out low, median, and high values of var2 (facet) and var3 (group)
  if(d(vars)==3){
    dat$facet <- str_extract(dat$facet, '0\\.[0-9]+')
    facetVals <- dat$facet %>% unique() %>% num() %>% sort()
    dat %<>% mutate(facet=num(facet))
    dat$facet[dat$facet==facetVals[1]] <- 'Low'
    dat$facet[dat$facet==facetVals[2]] <- 'Median'
    dat$facet[dat$facet==facetVals[3]] <- 'High'
    dat %<>% mutate(facet = factor(facet, levels=c('Low', 'Median', 'High')))
  }
  groupVals <- sort(num(unique(dat$group)))
  
  # Rename values of +/- 1 SD for plotting
  dat %<>% mutate(group=num(group))
  dat$group[dat$group==groupVals[1]] <- 'Low'
  dat$group[dat$group==groupVals[2]] <- 'Median'
  dat$group[dat$group==groupVals[3]] <- 'High'
  
  dat %<>%
    mutate(group = factor(group, levels=c('Low', 'Median', 'High'))) 
  
  # Map x-axis values (+/- 1 SD from the median if SD=TRUE, but defaults to min/max)
  if(sd){
    xVar <- str_replace(vars[1], ' \\[all\\]', '')
    xMed <- median(mdat[,xVar], na.rm=TRUE)
    xSD <- sd(mdat[,xVar], na.rm=TRUE)
    xLo <- xMed-xSD
    xHi <- xMed+xSD
    matchVal <- function(vec, val){which.min(abs(vec-val))}
    xVals <- sapply(c(xLo, xHi), function(x){ dat$x[matchVal(dat$x, x)] })
  } else {
    xVals <- c(min(dat$x), max(dat$x))
  }
  
  
  # Filter to +/- 1 SD from the median for the x axis
  dat <- filter(dat, x %in% xVals)
  
  # Compute changes from low to high values of x axis
  if(d(vars)==3){
    pdat <- expand.grid(group=unique(dat$group), facet=unique(dat$facet))
  } else {
    pdat <- expand.grid(group=unique(dat$group))
  }
  pdat %<>%
    mutate(delta=NA) %>%
    mutate(lo=NA) %>%
    mutate(hi=NA)
  for(ii in 1:nrow(pdat)){
    if(d(vars)==3){
      df <- filter(dat, group==pdat$group[ii] & facet==pdat$facet[ii]) # Filter to relevant combination of maj and min spread
    } else {
      df <- filter(dat, group==pdat$group[ii]) # Filter to relevant combination of maj and min spread
    }
    pdat$delta[ii] <- round(df$predicted[df$x==max(df$x)] - df$predicted[df$x==min(df$x)], 3) #Calculate the change
    pdat$lo[ii] <- round(df$conf.low[df$x==max(df$x)]-df$conf.high[df$x==min(df$x)], 3)
    pdat$hi[ii] <- round(df$conf.high[df$x==max(df$x)]-df$conf.low[df$x==min(df$x)], 3)
  }
  
  # Return
  return(pdat)
  
}

# Run
pAll <- predDat(model=m.twoint, vars=c('medDiff [all]', 'majSD', 'minSD'))
pMin <- predDat(model=m.min, vars=c('medDiff [all]', 'majSD', 'minSD'))

#Plot together
# Set up plot parameters
pdf("~/Dropbox/Projects/Minority Party/Graphics/probVoteChange_3way_passageDV.pdf", width=10, height=5)
par(mfrow=c(1,2))

# Plot all bills
breaks <- (nrow(pAll)/3)-1
yVals <- c(pAll$lo, pAll$hi) %>% 
  range() %>%
  sapply(., function(x){ ifelse(x < 0, round_any(x, 0.1, f=floor), round_any(x, 0.1, f=ceiling)) })
plot(1:(nrow(pAll)+2*breaks), 1:(nrow(pAll)+(2*breaks)), type='n', 
     ylab='Change in Probability of Vote', xlab='Majority Spread', ylim=yVals, xaxt='n', las=2, main='All Bills')
abline(h=0, lty=2)
for(ii in 1:nrow(pAll)){
  colIndex <- which(c('Low', 'Median', 'High')==pAll$group[ii])
  if(ii %in% 1:3){
    points(ii, pAll$delta[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pAll$lo[ii], pAll$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pAll$delta[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pAll$lo[ii], pAll$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pAll$delta[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pAll$lo[ii], pAll$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Minority Spread')

# Plot minority sponsored bills
breaks <- (nrow(pMin)/3)-1
yVals <- c(pMin$lo, pMin$hi) %>% 
  range() %>%
  sapply(., function(x){ ifelse(x < 0, round_any(x, 0.1, f=floor), round_any(x, 0.1, f=ceiling)) })
plot(1:(nrow(pMin)+2*breaks), 1:(nrow(pMin)+(2*breaks)), type='n', 
     ylab='Change in Probability of Vote', xlab='Majority Spread', ylim=yVals, xaxt='n', las=2, main='Minority-Sponsored Bills')
abline(h=0, lty=2)
for(ii in 1:nrow(pMin)){
  colIndex <- which(c('Low', 'Median', 'High')==pMin$group[ii])
  if(ii %in% 1:3){
    points(ii, pMin$delta[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pMin$lo[ii], pMin$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pMin$delta[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pMin$lo[ii], pMin$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pMin$delta[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pMin$lo[ii], pMin$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Minority Spread')

# Turn off graphics device
par(mfrow=c(1,1))
dev.off()


